function [ys,params,check] = markov2_perfect_steadystate(ys,exo,M_,options_)
% function [ys,params,check] = NK_baseline_steadystate(ys,exo,M_,options_)
% computes the steady state for the NK_baseline.mod and uses a numerical

% read out parameters to access them with their name
NumberOfParameters = M_.param_nbr;
for ii = 1:NumberOfParameters
  paramname = M_.param_names{ii};
  eval([ paramname ' = M_.params(' int2str(ii) ');']);
end
% initialize indicator
check = 0;



%%%%%%%%%% Don't change the above code %%%%%%%%%%%%
global theta dbar  L gL alpha rbest rho sigma ns n betas d alphas fspace_smolyak ilinear ilog_rL delta beta


% private eqm
r_private = rbest;
n  = ncountries;
ns = nsectors;
ilog_rL = 1; % use exp(rL)/(1+exp(rL))

dij = dbar;
d=dij*ones(n); %dij is (1+iceberg) cost producing in j selling to i
d=d-(dij-1)*eye(n);


for i=1:n
    eval(['L(',num2str(i),',1)=L',num2str(i),';'])
    for j=1:ns
        eval(['alpha(',num2str(i),',',num2str(j),')=alpha',num2str(i),num2str(j),';'])
       % eval(['sa(',num2str(i),',',num2str(j),')=sa',num2str(i),num2str(j),';'])
       % alpha(i,j) = alpha1(i,j)*exp(sa(i,j));
    end
end

for j=1:ns
    eval(['betas(1,',num2str(j),')=betas',num2str(j),';'])
end

%% solve for the steady state of gov assuming dG2/dT=0!!!
% This is not the true steady state, it's only for us to get initial guess
% --- steady state ---


options = optimoptions('fsolve', 'TolFun', 1e-12);
options.Display = 'off';

xguess(1:2) = 1.1;
xguess(3) = 0.2;

[xval,fval_sol,eflag_ss] = fsolve(@(xval) solve_1s_2c_gov_ss_nogL(xval),xguess,options); %iter


[Tss,xss,Pss,pi_ss, fval_ss] = eval_1s_2c_gov_ss_nogL(xval);

xval_ss= xval;
wss = xval(1:n);
extax_ss = xval(3);
taf_ss = 0;
rss = rbest;

% construct M
for i=1:n
    GTss(i)= log(xss(i)^(-sigma)*Pss(i)^(sigma-1)*wss(i)/alpha(i));
    Mss(i) = wss(i)/alpha(i);
end

uc = xss(1)^(-sigma);

for i=2:n
    for j=1:ns
        for ii=1:n
            for jj=1:ns
                eval(['g0=gcoef0_',num2str(i),num2str(j),num2str(ii),num2str(jj),';']);
                dGT(i,ii) = g0*exp(GTss(i))/Tss(ii);
            end
        end
    end
end
               
gammaT_ss = wss(2)/alpha(2)*(1+theta)/theta*pi_ss(2,1)*uc*extax_ss/(1+extax_ss);
mu = gammaT_ss*alpha(2);

%% write to our variable names in dynare
for j=1:ns
    eval(['wedge',num2str(j),'=0;'])
end

rL = rbest*L;
Lp_ss=(1-rbest)*L;


eval(['bterm= 0;'])
for i=1:n
    eval(['x',num2str(i),'=log(xss(',num2str(i),'));'])
    eval(['w',num2str(i),'=log(wss(',num2str(i),'));'])
    eval(['P',num2str(i),'=log(Pss(',num2str(i),'));'])
    eval(['riskfree',num2str(i),'=log(1/beta);']);
    if (i>1)
        eval(['mu',num2str(i),'=mu;'])
    end

    for j=1:ns
        eval(['T',num2str(i),num2str(j),'=log(Tss(',num2str(i),'));'])
        eval(['rr',num2str(i),num2str(j),'=log(rbest);'])
        eval(['Lr',num2str(i),num2str(j),'=log(rL(',num2str(i),'));'])
        eval(['Lp',num2str(i),num2str(j),'=log(Lp_ss(',num2str(i),'));'])
        eval(['M',num2str(i),num2str(j),'=Mss(',num2str(i),');'])
        eval(['GT',num2str(i),num2str(j),'=GTss(',num2str(i),');'])
        eval(['sa',num2str(i),num2str(j),'=0;'])

        if (i>1)
            eval(['gammar',num2str(i),num2str(j),'= 0;'])
            eval(['taux',num2str(i),num2str(j),'= extax_ss;'])
            eval(['taf',num2str(i),num2str(j),'= taf_ss;'])
            eval(['gammaT',num2str(i),num2str(j),'= gammaT_ss;'])
            
            
        end


        for ii=1:n
            eval(['pi',num2str(i),num2str(ii),num2str(j),'= pi_ss(',num2str(i),',',num2str(ii),');'])
        end
      

%             for ii=1:n
%                 for jj=1:ns
%                     eval(['dGT',num2str(i),num2str(j),num2str(ii),num2str(jj),'=gcoef0_',num2str(i),num2str(j),...
%                         '+gcoef1_',num2str(i),num2str(j),...
%                         '*Tss(',num2str(i),',',num2str(j),')+gcoef2_',num2str(i),num2str(j),'*Tss(',num2str(i),',',num2str(j),')^2;']);
%                 end
%             end
        if (i>1)
            for ii=1:n
                for jj=1:ns
                    
                    eval(['dGT',num2str(i),num2str(j),num2str(ii),num2str(jj),...
                           '=dGT(',num2str(i),',',num2str(ii),');']);
                  
                end
            end

        end

    end
end

%% end own model equations, below from dynare

params=NaN(NumberOfParameters,1);
for iter = 1:length(M_.params) %update parameters set in the file
  eval([ 'params(' num2str(iter) ') = ' M_.param_names{iter} ';' ])
end

NumberOfEndogenousVariables = M_.orig_endo_nbr; %auxiliary variables are set automatically
for ii = 1:NumberOfEndogenousVariables
  varname = M_.endo_names{ii};
  eval(['ys(' int2str(ii) ') = ' varname ';']);
end
